clear all, close all, clc
global lamda sigma nu betta taumax depr c1 r Fkk T c c0 kss css c1ss rss Fkkss yss kstart kistart rnext c1next Fkknext

N=2;
betta=0.96;
sigma=3;
eta=3;
z=2.5;
rho=0.95;
xi=1;
depr=1;
taumax=0.1;
n=1;
lamdaBSz=1-(betta*(1-taumax)*z^(1-sigma)*rho)^(1/sigma);

%SS
kss=sqrt((((1-lamdaBSz)*z)^2-rho)/(1-rho))
yss=z*(rho/kss^2+(1-rho))^(1/(1-sigma));
css=lamdaBSz*yss
rss=rho*z^(1-eta)*(yss/kss)^eta;
Fkkss=-eta*rho*(1-rho)*z^(2-2*eta)*yss^(2*eta-1)*kss^(-eta-1);   
nu=1-rss*taumax/(1-depr+rss);

T=200;
k=zeros(1,T);
kstart=1;
kistart=[0 1];
k(1)=(1-lamdaBSz)*z*(rho+(1-rho)*kstart^2)^(-1/2)*kstart;
for i=2:T
    k(i)=(1-lamdaBSz)*z*(rho+(1-rho)*k(i-1)^2)^(-1/2)*k(i-1);
end
y(1)=z*(rho/kstart^2+(1-rho))^(1/(1-sigma));
y(2:T)=z*(rho./k(1:T-1).^2+(1-rho)).^(1/(1-sigma));
c=lamdaBSz*y;
r(1)=rho*z^(1-eta)*(y(1)./kstart).^eta;
r(2:T)=rho*z^(1-eta)*(y(2:T)./k(1:T-1)).^eta;
Fkk(1)=-eta*rho*(1-rho)*z^(2-2*eta)*y(1)^(2*eta-1)*kstart^(-eta-1);
Fkk(2:T)=-eta*rho*(1-rho)*z^(2-2*eta)*y(2:T).^(2*eta-1).*k(1:T-1).^(-eta-1);
cnext = [c(2:T) css];
rnext = [r(2:T) rss];
c.^(-sigma)- cnext.^(-sigma)*betta.*(1+rnext*(1-taumax)-depr); %Euler
Fkknext = [Fkk(2:end) Fkkss];

kBsz=[kstart k(1:23)];

%check lamda
Dai=-1;
y0=z;
c0=lamdaBSz*y0;
A=(y0/c0)^(sigma)*(1-taumax)*rho*z^(1-sigma);
num=lamdaBSz^(1-sigma)*rho;
denom=1-betta*rho*((1-lamdaBSz)*z)^(1-sigma);
xx1=-0.08*num/denom;
part1=lamdaBSz^(1-sigma)*(1-rho);
part2=betta*xx1/(-0.08)*((1-lamdaBSz)*z)^(1-sigma)*(1-rho);
xx2=-0.08*(part1+part2)/(1-betta);
V=xx1+xx2;
lamda=1+A/V/(1-sigma)*Dai %alpha^i in BSz; using their method

options=optimset('TolFun',1e-20,'TolX',1e-20,'Disp','iter','MaxFunEvals',1e5);
lamda=fsolve('findLam',0.5,options)

c1ss=lamda*css;
c1=lamda*c;
c1next = [c1(2:T) c1ss];
c10=lamda*c0;

for i=1:T;
    X(i) = 1; %gamma
end;
Delta1=1;
deductible=1;
X(T+1:T+2)=[Delta1 deductible];

resid3v_taumaxIIIA(X);

options = optimset('MaxFunEvals',10000000,'MaxIter',200,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
[X,check] = fsolve('resid3v_taumaxIIIA',X,options)
sum(resid3v_taumaxIIIA(X).^2)

gamma=X(1:T);
Delta1 = X(T+1)
Delta2 = - Delta1;
deductible = X(T+2)

options=optimset('TolFun',1e-20,'TolX',1e-20,'Disp','off','MaxFunEvals',1e5); 
gammass=fsolve('findgammass',0.5,options);
muss=lamda*((1+(Delta1+Delta2/lamda)*(1-sigma))*c1ss^(-sigma)-sigma*c1ss^(-sigma-1)*gammass*(1-rss*(1-taumax)));
Omegabar=1+(Delta1+Delta2/lamda)*(1-sigma)

gammanext = [X(2:T) gammass];

tauK0=taumax;
tauKt(1:T) = taumax;
tauKt = [tauKt taumax];

%period 0
c0 = c1(1);
r0 = r(1);

%periods t>0
ct = c1(2:T);
rt = r(2:T);

transfer=0;

mu(1)=lamda*((1+(Delta1+Delta2/lamda)*(1-sigma))*c0^(-sigma)...
    -sigma*gamma(1)*c0^(-sigma-1)...
    + sigma*c0^(-sigma-1)*((Delta1*kistart(1)+Delta2*kistart(2))*(1+r0*(1-tauK0)-depr)+Delta1*(deductible-transfer)+Delta2*(deductible+transfer))); %FOC for consumption at time 0
mu(2:T)=lamda*((1+(Delta1+Delta2/lamda)*(1-sigma))*ct.^(-sigma)...
    +gamma(2:T)*(-sigma).*ct.^(-sigma-1)...
    -gamma(1:T-1)*(-sigma).*ct.^(-sigma-1).*(1+rt.*(1-taumax)-depr)); %FOC for consumption at t>0

figure
subplot(121)
plot(gamma,'-k','LineWidth',2)
ylim([-0.05 0.3])
yline(0,':k','LineWidth',1.3);
h1 = xlabel('Time'); % Create label 
pos = get(h1,'pos'); % Read position [x y z] 
set(h1,'pos',pos+[0 0 0]) % Move label
title('$\gamma$','Interpreter','latex','FontSize',16);

subplot(122)
plot(mu,'-k','LineWidth',2)
ylim([-0.05 0.3])
yline(0,':k','LineWidth',1.3);
h1 = xlabel('Time'); % Create label 
pos = get(h1,'pos'); % Read position [x y z] 
set(h1,'pos',pos+[0 0 0]) % Move label
title('$\mu$','Interpreter','latex','FontSize',16);
